/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;
using System.Collections.Generic;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.LinearAlgebra.Solvers.Preconditioners;
using dnAnalytics.Statistics.Distributions;
using dnAnalytics.Math;
using System.Diagnostics;

namespace dnAnalytics.LinearAlgebra.Solvers.Iterative
{
    /// <summary>
    /// A Multiple-Lanczos Bi-Conjugate Gradient stabilized iterative matrix solver.
    /// </summary>
    /// <remarks>
    /// <para>
    /// The Multiple-Lanczos Bi-Conjugate Gradient stabilized (ML(k)-BiCGStab) solver is an 'improvement'
    /// of the standard BiCgStab solver. 
    /// </para>
    /// <para>
    /// The algorithm was taken from: <br/>
    /// ML(k)BiCGSTAB: A BiCGSTAB variant based on multiple Lanczos starting vectors
    /// <br/>
    /// Man-chung Yeung and Tony F. Chan
    /// <br/>
    /// SIAM Journal of Scientific Computing
    /// <br/>
    /// Volume 21, Number 4, pp. 1263 - 1290
    /// </para>
    /// <para>
    /// The example code below provides an indication of the possible use of the
    /// solver.
    /// </para>
    /// </remarks>
    /// <include file='../../../../examples/LinearAlgebra/Solvers/Iterative/MlkBiCgStab.xml' path='example'/> 
    public sealed class MlkBiCgStab : IIterativeSolver
    {
        /// <summary>
        /// The default number of starting vectors.
        /// </summary>
        private const int sm_DefaultNumberOfStartingVectors = 50;
        
        /// <summary>
        /// The status used if there is no status, i.e. the solver hasn't run yet and there is no
        /// iterator.
        /// </summary>
        private static readonly ICalculationStatus sm_DefaultStatus = new CalculationIndetermined();

        /// <summary>
        /// The preconditioner that will be used. Can be set to <c>null</c>, in which case the default
        /// pre-conditioner will be used.
        /// </summary>
        private IPreConditioner m_Preconditioner;

        /// <summary>
        /// The iterative process controller.
        /// </summary>
        private IIterator m_Iterator;

        /// <summary>
        /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
        /// </summary>
        private IList<Vector> m_StartingVectors;

        /// <summary>
        /// The number of starting vectors used by the algorithm
        /// </summary>
        private int m_NumberOfStartingVectors = sm_DefaultNumberOfStartingVectors;

        /// <summary>
        /// Indicates if the user has stopped the solver.
        /// </summary>
        private bool m_HasBeenStopped;

        /// <summary>
        /// Initializes a new instance of the <c>MlkBiCgStab</c> class.
        /// </summary>
        /// <remarks>
        /// When using this constructor the solver will use the <see cref="IIterator"/> with
        /// the standard settings and a default preconditioner.
        /// </remarks>
        public MlkBiCgStab() : this(null, null)
        {}

        /// <summary>
        /// Initializes a new instance of the <c>MlkBiCgStab</c> class.
        /// </summary>
        /// <remarks>
        /// <para>
        /// When using this constructor the solver will use a default preconditioner.
        /// </para>
        /// <para>
        /// The main advantages of using a user defined <see cref="IIterator"/> are:
        /// <list type="number">
        /// <item>It is possible to set the desired convergence limits.</item>
        /// <item>
        /// It is possible to check the reason for which the solver finished 
        /// the iterative procedure by calling the <see cref="IIterator.Status"/> property.
        /// </item>
        /// </list>
        /// </para>
        /// </remarks>
        /// <param name="iterator">
        ///   The <see cref="IIterator"/> that will be used to monitor the iterative process.
        /// </param>
        public MlkBiCgStab(IIterator iterator) : this(null, iterator)
        {}

        /// <summary>
        /// Initializes a new instance of the <c>MlkBiCgStab</c> class.
        /// </summary>
        /// <remarks>
        /// When using this constructor the solver will use the <see cref="IIterator"/> with
        /// the standard settings.
        /// </remarks>
        /// <param name="preconditioner">
        /// The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.
        /// </param>
        public MlkBiCgStab(IPreConditioner preconditioner) : this(preconditioner, null)
        {}

        /// <summary>
        /// Initializes a new instance of the <c>MlkBiCgStab</c> class.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The main advantages of using a user defined <see cref="IIterator"/> are:
        /// <list type="number">
        /// <item>It is possible to set the desired convergence limits.</item>
        /// <item>
        /// It is possible to check the reason for which the solver finished 
        /// the iterative procedure by calling the <see cref="IIterator.Status"/> property.
        /// </item>
        /// </list>
        /// </para>
        /// </remarks>
        /// <param name="preconditioner">
        /// The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.
        /// </param>
        /// <param name="iterator">
        /// The <see cref="IIterator"/> that will be used to monitor the iterative process.
        /// </param>
        public MlkBiCgStab(IPreConditioner preconditioner, IIterator iterator)
        {
            m_Iterator = iterator;
            m_Preconditioner = preconditioner;
        }

        /// <summary>
        /// Gets or sets the number of starting vectors.
        /// </summary>
        /// <remarks>
        ///   Must be larger than 1 and smaller than the number of variables in the matrix that 
        ///   for which this solver will be used.
        /// </remarks>
        public int NumberOfStartingVectors
        {
            [DebuggerStepThrough]
            get
            {
                return m_NumberOfStartingVectors;
            }
            [DebuggerStepThrough]
            set
            {
                if (value <= 1)
                {
                    throw new ArgumentOutOfRangeException("value");
                }

                m_NumberOfStartingVectors = value;
            }
        }

        /// <summary>
        /// Resets the number of starting vectors to the default value.
        /// </summary>
        public void ResetNumberOfStartingVectors()
        {
            m_NumberOfStartingVectors = sm_DefaultNumberOfStartingVectors;
        }

        /// <summary>
        /// Sets the <c>IPreConditioner</c> that will be used to precondition the iterative process.
        /// </summary>
        /// <param name="preconditioner">The preconditioner.</param>
        public void SetPreconditioner(IPreConditioner preconditioner)
        {
            m_Preconditioner = preconditioner;
        }

        /// <summary>
        /// Sets the <c>IIterator</c> that will be used to track the iterative process.
        /// </summary>
        /// <param name="iterator">The iterator.</param>
        public void SetIterator(IIterator iterator)
        {
            m_Iterator = iterator;
        }

        /// <summary>
        /// Gets or sets a series of orthonormal vectors which will be used as basis for the 
        /// Krylov sub-space.
        /// </summary>
        /// <remarks>
        /// </remarks>
        public IList<Vector> StartingVectors
        {
            [DebuggerStepThrough]
            get
            {
                return m_StartingVectors;
            }
            [DebuggerStepThrough]
            set
            {
                if ((value == null) || (value.Count == 0))
                {
                    m_StartingVectors = null;
                }
                else
                {
                    m_StartingVectors = value;
                }
            }
        }

        /// <summary>
        /// Gets the status of the iteration once the calculation is finished.
        /// </summary>
        public ICalculationStatus IterationResult
        {
            [DebuggerStepThrough]
            get
            {
                return (m_Iterator != null) ? m_Iterator.Status : sm_DefaultStatus;
            }
        }

        /// <summary>
        /// Stops the solve process. 
        /// </summary>
        /// <remarks>
        /// Note that it may take an indetermined amount of time for the solver to actually
        /// stop the process.
        /// </remarks>
        public void StopSolve()
        {
            m_HasBeenStopped = true;
        }

        /// <summary>
        /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
        /// solution vector and x is the unknown vector.
        /// </summary>
        /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
        /// <param name="vector">The solution vector, <c>b</c>.</param>
        /// <returns>The result vector, <c>x</c>.</returns>
        public Vector Solve(Matrix matrix, Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException();
            }

            Vector result = new DenseVector(matrix.Rows);
            Solve(matrix, vector, result);
            return result;
        }

        /// <summary>
        /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
        /// solution vector and x is the unknown vector.
        /// </summary>
        /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
        /// <param name="input">The solution vector, <c>b</c></param>
        /// <param name="result">The result vector, <c>x</c></param>
        public void Solve(Matrix matrix, Vector input, Vector result)
        {
            // If we were stopped before, we are no longer
            // We're doing this at the start of the method to ensure
            // that we can use these fields immediately.
            {
                m_HasBeenStopped = false;
            }

            // Error checks
            {
                if (matrix == null)
                {
                    throw new ArgumentNullException("matrix");
                }

                if (matrix.Rows != matrix.Columns)
                {
                    throw new MatrixNotSquareException();
                }

                if (input == null)
                {
                    throw new ArgumentNullException("input");
                }

                if (result == null)
                {
                    throw new ArgumentNullException("result");
                }

                if (result.Count != input.Count)
                {
                    throw new NotConformableException();
                }
            }

            // Initialize the solver fields
            {
                // Set the convergence monitor
                if (m_Iterator == null)
                {
                    m_Iterator = Iterator.CreateDefault();
                }

                if (m_Preconditioner == null)
                {
                    m_Preconditioner = new UnitPreconditioner();
                }
                m_Preconditioner.Initialize(matrix);
            }


            // Choose an initial guess x_0
            // Take x_0 = 0
            Vector xTemp = new DenseVector(input.Count);

            // Choose k vectors q_1, q_2, ..., q_k
            {
                // Build a new set if:
                // a) the stored set doesn't exist (i.e. == null)
                // b) Is of an incorrect length (i.e. too long)
                // c) The vectors are of an incorrect length (i.e. too long or too short)
                bool useOld = false;
                if (m_StartingVectors != null)
                {
                    // We don't accept collections with zero starting vectors so ...
                    if (m_StartingVectors.Count <= NumberOfStartingVectorsToCreate(m_NumberOfStartingVectors, input.Count))
                    {
                        // Only check the first vector for sizing. If that matches we assume the 
                        // other vectors match too. If they don't the process will crash
                        if (m_StartingVectors[0].Count == input.Count)
                        {
                            useOld = true;
                        }
                    }
                }

                m_StartingVectors = (useOld) ? m_StartingVectors : CreateStartingVectors(m_NumberOfStartingVectors, input.Count);
            }

            // Store the number of starting vectors. Not really necessary but easier to type :)
            int k = m_StartingVectors.Count;

            // r_0 = b - Ax_0
            // This is basically a SAXPY so it could be made a lot faster
            Vector residuals = new DenseVector(matrix.Rows);
            CalculateTrueResidual(matrix, residuals, xTemp, input);

            // Define the temporary scalars
            double[] c = new double[k];
            double alpha = 0;
            double rho = 0;
            double beta = 0;

            // Define the temporary vectors
            Vector gTemp = new DenseVector(residuals.Count);

            Vector u = new DenseVector(residuals.Count);
            Vector uTemp = new DenseVector(residuals.Count);
            Vector temp = new DenseVector(residuals.Count);
            Vector temp1 = new DenseVector(residuals.Count);

            Vector zd = new DenseVector(residuals.Count);
            Vector zg = new DenseVector(residuals.Count);
            Vector zw = new DenseVector(residuals.Count);

            Vector[] d = CreateVectorArray(m_StartingVectors.Count, residuals.Count);
            // g_0 = r_0
            Vector[] g = CreateVectorArray(m_StartingVectors.Count, residuals.Count);
            residuals.CopyTo(g[k - 1]);

            Vector[] w = CreateVectorArray(m_StartingVectors.Count, residuals.Count);

            // FOR (j = 0, 1, 2 ....)
            int iterationNumber = 0;
            while (ShouldContinue(iterationNumber, xTemp, input, residuals))
            {
                // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k)
                m_Preconditioner.Approximate(g[k - 1], gTemp);

                // w_((j-1)k+k) = A g~_((j-1)k+k)
                matrix.Multiply(gTemp, w[k - 1]);

                // c_((j-1)k+k) = q^T_1 w_((j-1)k+k)
                c[k - 1] = m_StartingVectors[0].DotProduct(w[k - 1]);
                if (Precision.EqualsWithTolerance(c[k - 1], 0, Constants.UlpsForComparison))
                {
                    throw new IterativeSolverBrakedownException();
                }

                // alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k)
                alpha = m_StartingVectors[0].DotProduct(residuals) / c[k - 1];

                // u_(jk+1) = r_((j-1)k+k) - alpha_(jk+1) w_((j-1)k+k)
                residuals.AddScaledVector(-alpha, w[k - 1], u);

                // SOLVE M u~_(jk+1) = u_(jk+1)
                m_Preconditioner.Approximate(u, temp1);
                temp1.CopyTo(uTemp);

                // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2
                matrix.Multiply(temp1, temp);

                rho = temp.DotProduct();
                // If rho is zero then temp is a zero vector and we're probably
                // about to have zero residuals (i.e. an exact solution).
                // So set rho to 1.0 because in the next step it will turn to zero.
                if (Precision.EqualsWithTolerance(rho, 0, Constants.UlpsForComparison))
                {
                    rho = 1.0;
                }
                rho = -u.DotProduct(temp) / rho;

                // x_(jk+1) = x_((j-1)k_k) - rho_(j+1) u~_(jk+1) + alpha_(jk+1) g~_((j-1)k+k)
                xTemp.AddScaledVector(-rho, uTemp);
                xTemp.AddScaledVector(alpha, gTemp);

                // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1)
                u.CopyTo(residuals);

                // Reuse temp
                temp.Multiply(rho);
                residuals.Add(temp);

                // Check convergence and stop if we are converged.
                if (!ShouldContinue(iterationNumber, xTemp, input, residuals))
                {
                    // Calculate the true residual
                    CalculateTrueResidual(matrix, residuals, xTemp, input);

                    // Now recheck the convergence
                    if (!ShouldContinue(iterationNumber, xTemp, input, residuals))
                    {
                        // We're all good now.
                        // Exit from the while loop.
                        break;
                    }
                }

                // FOR (i = 1,2, ...., k)
                for (int i = 0; i < k; i++)
                {
                    // z_d = u_(jk+1)
                    u.CopyTo(zd);

                    // z_g = r_(jk+i)
                    residuals.CopyTo(zg);

                    // z_w = 0
                    zw.Clear();

                    // FOR (s = i, ...., k-1) AND j >= 1
                    if (iterationNumber >= 1)
                    {
                        for (int s = i; s < k - 1; s++)
                        {
                            // beta^(jk+i)_((j-1)k+s) = -q^t_(s+1) z_d / c_((j-1)k+s)
                            beta = -m_StartingVectors[s + 1].DotProduct(zd) / c[s];

                            // z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s)
                            zd.AddScaledVector(beta, d[s]);

                            // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s)
                            zg.AddScaledVector(beta, g[s]);

                            // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s)
                            zw.AddScaledVector(beta, w[s]);
                        }
                    }

                    beta = rho * c[k - 1];
                    if (Precision.EqualsWithTolerance(beta, 0, Constants.UlpsForComparison))
                    {
                        throw new IterativeSolverBrakedownException();
                    }

                    // beta^(jk+i)_((j-1)k+k) = -(q^T_1 (r_(jk+1) + rho_(j+1) z_w)) / (rho_(j+1) c_((j-1)k+k))
                    residuals.AddScaledVector(rho, zw, temp);
                    beta = -m_StartingVectors[0].DotProduct(temp) / (beta);

                    // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k)
                    zg.AddScaledVector(beta, g[k - 1]);

                    // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k))
                    zw.AddScaledVector(beta, w[k - 1]);
                    zw.Multiply(rho);

                    // z_d = r_(jk+i) + z_w
                    residuals.Add(zw, zd);

                    // FOR (s = 1, ... i - 1)
                    for (int s = 0; s < i - 1; s++)
                    {
                        // beta^(jk+i)_(jk+s) = -q^T_s+1 z_d / c_(jk+s)
                        beta = -m_StartingVectors[s + 1].DotProduct(zd) / c[s];

                        // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s)
                        zd.AddScaledVector(beta, d[s]);

                        // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s)
                        zg.AddScaledVector(beta, g[s]);
                    }

                    // d_(jk+i) = z_d - u_(jk+i)
                    zd.Subtract(u, d[i]);

                    // g_(jk+i) = z_g + z_w
                    zg.Add(zw, g[i]);

                    // IF (i < k - 1)
                    if (i < k - 1)
                    {
                        // c_(jk+1) = q^T_i+1 d_(jk+i)
                        c[i] = m_StartingVectors[i + 1].DotProduct(d[i]);
                        if (Precision.EqualsWithTolerance(c[i], 0, Constants.UlpsForComparison))
                        {
                            throw new IterativeSolverBrakedownException();
                        }

                        // alpha_(jk+i+1) = q^T_(i+1) u_(jk+i) / c_(jk+i)
                        alpha = m_StartingVectors[i + 1].DotProduct(u) / c[i];

                        // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i)
                        u.AddScaledVector(-alpha, d[i]);

                        // SOLVE M g~_(jk+i) = g_(jk+i)
                        m_Preconditioner.Approximate(g[i], gTemp);

                        // x_(jk+i+1) = x_(jk+i) + rho_(j+1) alpha_(jk+i+1) g~_(jk+i)
                        xTemp.AddScaledVector(rho * alpha, gTemp);

                        // w_(jk+i) = A g~_(jk+i)
                        matrix.Multiply(gTemp, w[i]);

                        // r_(jk+i+1) = r_(jk+i) - rho_(j+1) alpha_(jk+i+1) w_(jk+i)
                        residuals.AddScaledVector(-rho * alpha, w[i]);

                        // We can check the residuals here if they're close
                        if (!ShouldContinue(iterationNumber, xTemp, input, residuals))
                        {
                            // Recalculate the residuals and go round again. This is done to ensure that
                            // we have the proper residuals.
                            CalculateTrueResidual(matrix, residuals, xTemp, input);
                        }
                    }
                } // END ITERATION OVER i

                iterationNumber++;
            }

            // copy the temporary result to the real result vector
            xTemp.CopyTo(result);
        }

        private static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables)
        {
            // Create no more starting vectors than the size of the problem - 1
            return System.Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1));
        }

        /// <summary>
        /// Returns an array of starting vectors.
        /// </summary>
        /// <param name="maximumNumberOfStartingVectors">
        ///   The maximum number of starting vectors that should be created.
        /// </param>
        /// <param name="numberOfVariables">
        ///   The number of variables.
        /// </param>
        /// <returns>
        ///   An array with starting vectors. The array will never be larger than the
        ///   <paramref name="maximumNumberOfStartingVectors"/> but it may be smaller if
        ///   the <paramref name="numberOfVariables"/> is smaller than 
        ///   the <paramref name="maximumNumberOfStartingVectors"/>.
        /// </returns>
        private IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
        {
            // Create no more starting vectors than the size of the problem - 1
            // Get random values and then orthogonalize them with
            // modified Gramm - Schmidt
            int count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables);

            // Get a random set of samples based on the standard normal distribution with
            // mean = 0 and sd = 1
            Normal distribution = new Normal();
            distribution.RandomNumberGenerator = new System.Random();

            Matrix matrix = new DenseMatrix(numberOfVariables, count);
            for (int i = 0; i < matrix.Columns; i++)
            {
                double[] samples = distribution.Sample(matrix.Rows);
                
                // Set the column
                matrix.SetColumn(i, samples);
            }

            // Compute the orthogonalization.
            GramSchmidt gs = new GramSchmidt(matrix);
            Matrix orthogonalMatrix = gs.Q();

            // Now transfer this to vectors
            List<Vector> result = new List<Vector>();
            for (int i = 0; i < orthogonalMatrix.Columns; i++)
            {
                result.Add(orthogonalMatrix.GetColumn(i));
                // Normalize the result vector
                result[i].Multiply(1 / result[i].PNorm(2));
            }

            return result;
        }

        private Vector[] CreateVectorArray(int arraySize, int vectorSize)
        {
            Vector[] result = new Vector[arraySize];
            for (int i = 0; i < result.Length; i++)
            {
                result[i] = new DenseVector(vectorSize);
            }

            return result;
        }

        // Calculates the true residual of the matrix equation Ax = b according to:
        // residual = b - Ax
        private void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
        {
            // -Ax = residual
            matrix.Multiply(x, residual);
            residual.Negate();

            // residual + b
            residual.Add(b);
        }

        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
        {
            if (m_HasBeenStopped)
            {
                m_Iterator.IterationCancelled();
                return true;
            }

            m_Iterator.DetermineStatus(iterationNumber, result, source, residuals);
            ICalculationStatus status = m_Iterator.Status;

            // We stop if either:
            // - the user has stopped the calculation
            // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.)
            return (!status.TerminatesCalculation) && (!m_HasBeenStopped);
        }

        /// <summary>
        /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
        /// solution matrix and X is the unknown matrix.
        /// </summary>
        /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
        /// <param name="input">The solution matrix, <c>B</c>.</param>
        /// <returns>The result matrix, <c>X</c>.</returns>
        public Matrix Solve(Matrix matrix, Matrix input)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            Matrix result = matrix.CreateMatrix(input.Rows, input.Columns);
            Solve(matrix, input, result);
            return result;
        }

        /// <summary>
        /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
        /// solution matrix and X is the unknown matrix.
        /// </summary>
        /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
        /// <param name="input">The solution matrix, <c>B</c>.</param>
        /// <param name="result">The result matrix, <c>X</c></param>
        public void Solve(Matrix matrix, Matrix input, Matrix result)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (matrix.Rows != input.Rows || input.Rows != result.Rows || input.Columns != result.Columns)
            {
                throw new NotConformableException();
            }

            for (int column = 0; column < input.Columns; column++)
            {
                Vector solution = Solve(matrix, input.GetColumn(column));
                foreach (KeyValuePair<int, double> element in solution.GetIndexedEnumerator())
                {
                    result.ValueAt(element.Key, column, element.Value);
                }
            }
        }
    }
}
